Poravnanja - 2. deo¶
Algoritmi za izračunavanje poravnanja sekvenci, predstavljeni u prošlom delu, mogu se dodatno unaprediti. Jedno od unapređenja ogleda se u boljoj vremenskoj i prostornoj složenosti algoritama. Svi, do sada, predstavljeni algoritmi imaju kvadratnu složenost $O(nm)$, gde je $n$ dužina prve a $m$ dužina druge sekvence.
Jedan pokušaj optimizacije može biti smanjenje prostorne složenosti, pamćenjem samo rezultata iz prethodne iteracije. Ukoliko su tabele dinamičkog programiranja bile popunjavane po redovima, dovoljno je pamtiti tekući i prethodni red tabele, u suprotnom, pamtiće se tekuća i prethodna kolona. Ovim pristupom, na prvi pogled, gubi se mogućnost rekonstrukcije poravnanja - računa se samo konačni skor poravnanja. Vremenska složenost algoritma i dalje ostaje kvadratna.
Pametnim pristupom izračunavanju, moguće je modifikovati algoritme tako da njihova prostorna složenost postane linearna, $O(min\{m,n\})$ (vremenska i dalje kvadratna) a da pri tom rekonstrukcija poravnanja bude moguća.
Primer efikasnijeg pristupa biće prikazan pomoću Hirshberg-ovog algoritma za optimizaciju Needleman-Wunsch algoritma
Dodatni materijali: Hirshberg's algorithm
Za početak, potrebno je implementirati Needleman-Wunsch algoritam tako da skor računa u linearnoj prostornoj složenosti (red po red / kolona po kolonu), pamteći samo vrednosti iz prethodne i tekuće iteracije. Povratna vrednost algoritma biće poslednji red tabele dinamičkog programiranja.
Primer¶
Neka je potrebno izračunati skor globalnog poravnanja sekvenci V = AGTACGCA i W = TATGC
Tabela koju konstruiše Needleman-Wunsch algoritam je:
T A T G C
0 -2 -4 -6 -8 -10
A -2 -1 0 -2 -4 -6
G -4 -3 -2 -1 0 -2
T -6 -2 -4 0 -2 -1
A -8 -4 0 -2 -1 -3
C -10 -6 -2 -1 -3 1
G -12 -8 -4 -3 1 -1
C -14 -10 -6 -5 -1 3
A -16 -12 -8 -7 -3 1Dok je rezultat poslednjeg reda tabele, dobijen kao rezultat prethodno navedenog algoritma linearne prostorne složenosti:
[-16, -12, -8, -7, -3, 1]
Postupak na koji Hirshberg-ov algoritam omogućava rekonstrukciju poravnanja, u linearnoj prostornoj složenosti, zasniva se na podeli-pa-vladaj strategiji konstrukcije algoritama.
Za početak, umesto izračunavanja cele tabele, računa se poslednji red (ili kolona, u zavisnosti od pristupa implementaciji) polovine tabele dinamičkog programiranja.
Primer
AGTACGCA = AGTA + CGCA
T A T G C
0 -2 -4 -6 -8 -10
A -2 -1 0 -2 -4 -6
G -4 -3 -2 -1 0 -2
T -6 -2 -4 0 -2 -1
A -8 -4 0 -2 -1 -3Nakon izračunavanja poslednjeg reda prve (gornje) polovine tabele, rezultat je:
[-8, -4, 0, -2, -1, -3]Druga polovina tabele se takođe može rekonstruisati u linearnoj složenosti krećući se u suprotnom smeru (od donjeg desnog ugla ka gornjem levom). Kako bi se to ostvarilo, potrebno je izračunati poravnanja obrnutih sekvenci.
Primer
AGTACGCA = AGTA + CGCAUmesto poravnanja CGCA sa sekvencom TATGC, poravnaće se obrnuta sekvenca ACGC sa obrnutom sekvencom CGTAT. Rezultat ovog poravnanja je tabela:
C G T A T
0 -2 -4 -6 -8 -10
A -2 -1 -3 -5 -4 -6
C -4 0 -2 -4 -6 -5
G -6 -2 2 0 -2 -4
C -8 -4 0 1 -1 -3Poslednji red ove tabele, obrnut kako bi se ispravilo prethodno obrtanje sekvenci za potrebe izračunavanja, je:
[-8, -4, 0, 1, -1, -3] -> [-3, -1, 1, 0, -4, -8]Sada je potrebno sagledati dobijene rezultate iz prethodnih koraka. Poslednji redovi prve i druge polovine izračunavanja su:
[-8, -4, 0, -2, -1, -3]
[-3, -1, 1, 0, -4, -8]Sumiranjem izračunatih redova dobija se:
[-8, -4, 0, -2, -1, -3]
[-3, -1, 1, 0, -4, -8] +
--------------------------------
[-11, -5, 1, -2, -5, -11]Maksimum rezultujućeg niza predstavlja poziciju koja maksimizira poravnanja obe polovine, pa će i globalno poravnanje polaznih sekvenci sigurno prolaziti kroz tu poziciju.
[-11, -5, *1, -2, -5, -11]Neka je pozicija maksimuma označena indeksom j. Poravnanje se rekurzivno nastavlja na nivou polovina $V[0:\frac{n}{2}]$ sa $W[0:j]$ i $V[\frac{n}{2}:]$ sa $W[j:]$
Dolaskom do slučaja kada je deljenjem jedna od sekvenci postala prazna, izvršava se bazni slučaj rekurzije, koji podrazumeva da se preostali karakteri neprazne sekvence poravnavaju sa prazninom. U slučaju kada jedna od niski ima dužinu $1$, primenjuje se originalni Needleman-Wunsch algoritam.
Konačni rezultat poravnanja dobija se spajanjem rezultata iz listova stabla rekurzije.
Poravnanje sa afinim računanjem cene praznine¶
Umesto korišćenja jedne tabele dinamičkog poravnanja, za potrebe afinog poravnanja koristiće se tri tabele. Intuitivno, svaka od tabela odgovara kretanjem određenom operacijom (insercija, delecija, uparivanje karaktera). "Donja" tabela odgovara poravnavanju dodavanjem praznina u prvoj sekvenci, "srednja" tabela odgovara uparivanju karaktera i "gornja" poravnavanju dodavanjem praznina u drugoj sekvenci. Kroz "donju" tabelu dozvoljeno je kretanje (prenos rezultata) samo sa polja iznad (insercija), kroz "srednju" tabelu dozvoljeno je kretanje sa prethodne dijagonalne pozicije dok se kroz "gornju" tabelu može kretati samo sa leve strane (delecija). Skok sa srednje tabele na bilo donju ili gornju nosi cenu $\sigma$ dok svaki korak kroz gornju ili donju tabelu košta dodatno $\epsilon$. Cena kretanje kroz srednju tabelu zavisi od uparivanja odgovarajućih karaktera, ali se razmatra i moguć dolazak na srednju tabelu sa neke od preostale dve. U svakom koraku bira se rezultat koji maksimizuje rezultat.
Prvi redovi i kolone tabela inicijalizovani su na nulu
Implementacije algoritama¶
class Alignments:
def __init__(self,
gap_penalty,
match_score,
missmatch_score,
aff_gap_start_penalty,
aff_next_gap_penalty,
nucleotide_match_scores = None,
aa_match_scores = None
):
self.gap_penalty = gap_penalty
self.match_score = match_score
self.missmatch_score = missmatch_score
self.nucleotide_match_scores = nucleotide_match_scores
self.aa_match_scores = aa_match_scores
self.sigma = aff_gap_start_penalty
self.eps = aff_next_gap_penalty
# Metod izračunava skor poravnanja dva karaktera
def score(self, c1, c2, c_type):
if c_type == 'nucleotide' and self.nucleotide_match_scores != None:
return nucleotide_match_scores[c1][c2]
if c_type == 'aminoacid' and self.aa_match_scores != None:
return nucleotide_match_scores[c1][c2]
if c1 == c2:
return self.match_score
else:
return self.missmatch_score
# Metod rekonstruiše poravnanje
def backtracking(self,
backtrack,
s1,
s2,
D = None,
stop_value = None,
start_i = None,
start_j = None
):
if start_i != None:
i = start_i
else:
i = len(s1)
if start_j != None:
j = start_j
else:
j = len(s2)
s1_alignment = ''
s2_alignment = ''
while backtrack[i][j] != None and (D == None or stop_value == None or (D[i][j] != stop_value)):
(next_i, next_j) = backtrack[i][j]
if next_i == i - 1 and next_j == j:
s1_alignment = s1[i - 1] + s1_alignment
s2_alignment = '-' + s2_alignment
elif next_i == i - 1 and next_j == j - 1:
s1_alignment = s1[i - 1] + s1_alignment
s2_alignment = s2[j - 1] + s2_alignment
else:
s1_alignment = '-' + s1_alignment
s2_alignment = s2[j - 1] + s2_alignment
i = next_i
j = next_j
return s1_alignment, s2_alignment
# Metod izračunava edit rastojanje i odgovarajuće poravnanje sekvenci s1 i s2
def edit(self, s1, s2):
n = len(s1)
m = len(s2)
# Inicijalizacija
D = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack[0][0] = None
# Bazni slučajevi
for i in range(1, n + 1):
D[i][0] = i
backtrack[i][0] = (i - 1, 0)
for j in range(1, m + 1):
D[0][j] = j
backtrack[0][j] = (0, j - 1)
# Rekurentne veze
for i in range(1, n + 1):
for j in range(1, m + 1):
from_top = D[i - 1][j] + 1
from_left = D[i][j - 1] + 1
from_diagonal = D[i - 1][j - 1] + int(s1[i - 1] != s2[j - 1])
D[i][j] = min(from_top, from_left, from_diagonal)
if D[i][j] == from_diagonal:
backtrack[i][j] = (i - 1, j - 1)
elif D[i][j] == from_top:
backtrack[i][j] = (i - 1, j)
else:
backtrack[i][j] = (i, j - 1)
# Rekonstrukcija poravnanja
s1_alignment, s2_alignment = self.backtracking(backtrack, s1, s2)
return s1_alignment, s2_alignment, D[n][m]
# Metod izračunava globalno poravnanje sekvenci s1 i s2
def needleman_wunsch(self, s1, s2, c_type = 'text'):
n = len(s1)
m = len(s2)
# Inicijalizacija
D = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack[0][0] = None
# Bazni slučajevi
for i in range(1, n + 1):
D[i][0] = i * self.gap_penalty
backtrack[i][0] = (i - 1, 0)
for j in range(1, m + 1):
D[0][j] = j * self.gap_penalty
backtrack[0][j] = (0, j - 1)
# Rekurentne veze
for i in range(1, n + 1):
for j in range(1, m + 1):
from_top = D[i - 1][j] + self.gap_penalty
from_left = D[i][j - 1] + self.gap_penalty
from_diagonal = D[i - 1][j - 1] + self.score(s1[i - 1], s2[j - 1], c_type)
D[i][j] = max(from_top, from_left, from_diagonal)
if D[i][j] == from_diagonal:
backtrack[i][j] = (i - 1, j - 1)
elif D[i][j] == from_top:
backtrack[i][j] = (i - 1, j)
else:
backtrack[i][j] = (i, j - 1)
# Rekonstrukcija poravnanja
s1_alignment, s2_alignment = self.backtracking(backtrack, s1, s2)
return s1_alignment, s2_alignment, D[n][m]
# Algoritam izračunava lokalno poravnanje sekvenci s1 i s2
def smith_waterman(self, s1, s2, c_type = 'text'):
n = len(s1)
m = len(s2)
# Inicijalizacija
D = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack[0][0] = None
# Bazni slučajevi
for i in range(1, n + 1):
D[i][0] = 0
backtrack[i][0] = (i - 1, 0)
for j in range(1, m + 1):
D[0][j] = 0
backtrack[0][j] = (0, j - 1)
max_i = 0
max_j = 0
max_score = 0
# Rekurentne veze
for i in range(1, n + 1):
for j in range(1, m + 1):
from_top = D[i - 1][j] + self.gap_penalty
from_left = D[i][j - 1] + self.gap_penalty
from_diagonal = D[i - 1][j - 1] + self.score(s1[i - 1], s2[j - 1], c_type)
D[i][j] = max(from_top, from_left, from_diagonal, 0)
if D[i][j] > max_score:
max_i = i
max_j = j
max_score = D[i][j]
if D[i][j] == from_diagonal:
backtrack[i][j] = (i - 1, j - 1)
elif D[i][j] == from_top:
backtrack[i][j] = (i - 1, j)
else:
backtrack[i][j] = (i, j - 1)
# Rekonstrukcija poravnanja
s1_alignment, s2_alignment = self.backtracking(backtrack, s1, s2, D, 0, max_i, max_j)
return s1_alignment, s2_alignment, D[n][m]
# Metod izračunava skor globalnog poravnanja sekvenci s1 i s2 u linearnom prostoru,
# kao povratnu vrednost vraća poslednji red tabele dinamičkog programiranja
def needleman_wunsch_last_line(self, s1, s2, c_type='text'):
n = len(s1)
m = len(s2)
D = [[0 for _ in range(m + 1)] for _ in range(2)]
for j in range(1, m + 1):
D[0][j] = j * self.gap_penalty
for i in range(1, n + 1):
D[1][0] = i * self.gap_penalty
for j in range(1, m + 1):
from_top = D[0][j] + self.gap_penalty
from_left = D[1][j - 1] + self.gap_penalty
from_diagonal = D[0][j - 1] + self.score(s1[i - 1], s2[j - 1], c_type)
D[1][j] = max(from_top, from_left, from_diagonal)
D[0][:] = D[1][:]
return D[1]
# Metod izračunava skor globalnog poravnanja sekvenci s1 i s2 u linearnom prostoru,
# uz rekonstrukciju poravnanja
def hirschberg(self, s1, s2, c_type='text'):
s1_alignment = ''
s2_alignment = ''
n = len(s1)
m = len(s2)
if n == 0:
s1_alignment = m * '-'
s2_alignment = s2
elif m == 0:
s1_alignment = s1
s2_alignment = n * '-'
if n == 1 or m == 1:
s1_alignment, s2_alignment, _ = self.needleman_wunsch(s1, s2, c_type)
else:
mid = n // 2
score_l = self.needleman_wunsch_last_line(s1[:mid], s2, c_type)
score_r = self.needleman_wunsch_last_line(s1[mid:][::-1], s2[::-1], c_type)
score_r.reverse()
max_score = float('-inf')
max_score_j = None
for j in range(m):
curr_score = score_l[j] + score_r[j]
if curr_score > max_score:
max_score = curr_score
max_score_j = j
s1_l_alignment, s2_l_alignment = self.hirschberg(s1[:mid], s2[:max_score_j], c_type)
s1_r_alignment, s2_r_alignment = self.hirschberg(s1[mid:], s2[max_score_j:], c_type)
s1_alignment = s1_l_alignment + s1_r_alignment
s2_alignment = s2_l_alignment + s2_r_alignment
return s1_alignment, s2_alignment
# Metod izračunava globalno poravnanje sekvenci s1 i s2 uz korišćenje
# afine funkcije za računanje kazne poravnanja sa prazninama
def affine_gap_score_alignment(self, s1, s2, c_type='text'):
n = len(s1)
m = len(s2)
D_lower = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
D_middle = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
D_upper = [[0 for _ in range(m + 1)] for _ in range(n + 1)]
backtrack = [[None for _ in range(m + 1)] for _ in range(n + 1)]
for i in range(1, n):
backtrack[i][0] = (i - 1, 0)
for j in range(1, m):
backtrack[0][j] = (0, j - 1)
for i in range(1, n + 1):
for j in range(1, m + 1):
D_lower[i][j] = max(
D_lower[i - 1][j] + self.eps,
D_middle[i - 1][j] + self.sigma
)
D_upper[i][j] = max(
D_upper[i][j - 1] + self.eps,
D_middle[i][j - 1] + self.sigma
)
D_middle[i][j] = max(
D_lower[i][j],
D_upper[i][j],
D_middle[i - 1][j - 1] + self.score(s1[i - 1], s2[j - 1], c_type)
)
if D_middle[i][j] == D_middle[i - 1][j - 1] + self.score(s1[i - 1], s2[j - 1], c_type):
backtrack[i][j] = (i - 1, j - 1)
elif D_middle[i][j] == D_upper[i][j]:
backtrack[i][j] = (i, j - 1)
else:
backtrack[i][j] = (i - 1, j)
s1_alignment, s2_alignment = self.backtracking(backtrack, s1, s2)
return s1_alignment, s2_alignment, D_middle[n][m]
nucleotide_match_scores = {
'A': {
'A': 1,
'T': -0.1,
'G': -0.3,
'C': -0.9
},
'T': {
'A': -0.1,
'T': 1,
'G': -0.3,
'C': -0.2
},
'C': {
'A': -0.9,
'T': -0.2,
'G': -0.3,
'C': 2
},
'G': {
'A': -0.3,
'T': -0.3,
'G': 1,
'C': -0.3
}
}
al = Alignments(
gap_penalty = -2,
match_score = 2,
missmatch_score = -1,
nucleotide_match_scores = nucleotide_match_scores,
aff_gap_start_penalty = -2,
aff_next_gap_penalty = -1
)
s1 = 'AGGCTAG'
s2 = 'CGAAGC'
s1_a, s2_a = al.hirschberg(s1, s2, 'text')
print(s1_a)
print(s2_a)
s1_a1, s2_a1, _ = al.affine_gap_score_alignment(s1, s2)
print(s1_a1)
print(s2_a1)